%
% blocksmooth.m -- perform block GS lex smoothing with overlapping blocks
%
%   input: u  -- current solution  
%          A  -- matrix (Nx*Ny) x (Nx*Ny)
%          f  -- rhs
%          Nx -- number of points in the x-direction
%          Ny -- number of points in the y-direction
%          bx -- block size in x direction
%          by -- block size in y direciton
%          sx -- skip size in x-direction, overlap size is bx-sx
%          sy -- skip size in y-direction  overlap size is by-sy
%          ns -- number of smoothing steps to apply
%
%   output: v -- smoothed solution
%
function v = blocksmooth(u,A,f,ns,Nx,Ny,bx,by,sx,sy)
  
  
  % initialize
  %
  v = u;
  
  % apply some krylov iterations before smoothing
  %
  %v = pcg(A,f,1e-10,2,[],[],u);
  %v = gmres(A,f,2,1e-10,1,[],[],u);
  %v = bicgstab(A,f,1e-10,2,[],[],u);

  for n=1:ns
    
    for j1=1:sy:Ny
      for i1=1:sx:Nx
            
        % form the restriction operator
        %
        i2 = min(i1 + (bx-1), Nx );
        j2 = min(j1 + (by-1), Ny );
        
        bt = (i2-i1+1)*(j2-j1+1);  % total block size
        I = zeros(bt,1);
        J = I;
        k = 1;
        for i3=i1:i2
          for j3=j1:j2
            I(k) = i3;
            J(k) = j3;
            k=k+1;
          end
        end
        
        K = sub2ind([Nx,Ny],J,I); % double check whether I or J 
                                  % should come first    
        R = spalloc(bt,Nx*Ny,bt);
        R(1:bt,K) = speye(bt);
        
        % compute the residual and restrict
        %
        RA = R*A;
        rr = R*f - RA*v;
        
           
        % form the restricted matrix, and compute the update
        %
        Ak = RA*R';
        du = Ak\rr;
                
        % update the solution
        %
        v = v + R'*du;
        
      end
    end
  end
  
